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Abstract. Massive stars drive strong winds that impact the surrounding interstellar medium, 
producing parsec-scale bubbles for isolated stars and superbubbles around young clusters. These 
bubbles can be observed across the electromagnetic spectrum, both the wind itself and the 
swept up interstellar gas. Runaway massive stars produce bow shocks that strongly compresses 
interstellar gas, producing bright infrared, optical and radio nebulae. With the detection of 
non-thermal radio emission from bow shocks, particle acceleration can now also be investigated. 
I review research on wind bubbles and bow shocks around massive stars, highlighting recent 
advances in infrared, radio and X-ray observations, and progress in multidimensional simula- 
tions of these nebulae. These advances enable quantitative comparisons between theory and 
observations and allow to test the importance of some physical processes such as thermal con- 
duction and Kelvin-Helmholtz instability in shaping nebulae and in constraining the energetics 
of stellar-wind feedback to the interstellar medium. 
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1. Models for spherical wind bubbles 


Massive stars have a strong effect on their surroundings through their intense radia- 
tion (especially extreme-UV, ionizing radiation), strong winds, eruptive explosions and 
supernova explosions at the end of their lives. Quantifying these feedback effects is im- 
portant for understanding the dynamical and chemical evolution of galaxies, and also 
the structure and dynamics of the interstellar medium (ISM) in our own Galaxy. On 
smaller scales, modelling of circumstellar nebulae can also give clues as to the evolution- 
ary history of some nearby massive stars. Comparison of models with observations can 
give constraints on different physical processes in astrophysical plasmas, such as particle 
acceleration and thermal conduction. 

In one of the first papers studying the effects of stellar winds on the ISM, 
(1966) proposed that the optical cavity in the Rosette Nebula around NGC2244 is main- 
tained by dynamical pressure of the strong and high-velocity winds of the massive stars 


in the central cluster. [Dyson & de Vries| (1972) built on this work, developing a theory 


for the dynamics of wind-blown bubbles. This was further generalised by Castor et al. 
(1975), and again in the classic paper by [Weaver et all (1977). The basic picture is that 
of a spherical wind expanding from a star or group of stars and displacing the interstellar 
gas. From the contact discontinuity between the two media, a reverse (or termination) 
shock propagates backwards towards the star generating a hot bubble of shocked coronal 
gas with temperature T ~ 10° — 108 K. Similarly a forward shock propagates into the 
undisturbed ISM as long as the bubble expands supersonically. The ISM is photoionized 
by extreme-UV radiation from the hot star(s), producing an H II region around the star 
that usually extends well beyond the wind bubble, but which may be trapped by the 


shocked ISM [Freyer et al 2003). 
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(1977) introduced thermal conduction, which smoothes out the contact 


discontinuity and produces significant quantities of gas in the temperature range from 
10°—10" K, that emits UV and thermal X-ray radiation. This can reproduce the observed 
[O v1] line emission that is observed, but tends to overpredict thermal X-rays (see, e.g., 
[Toalá et al.2016). The dynamical model of Weaver has the bubble radius, r, expanding 
with time, t, as r œ t3/5. This is the solution in the adiabatic limit, i.e., in the limit 
that the thermal energy in the wind bubble is not lost and can power the pressure-driven 
expansion of the bubble. In the momentum conserving limit, where the wind energy is 
efficiently radiated away, a simple dimensional analysis shows that r œ t!/2. If thermal 
conduction is efficient at transporting thermal energy out of the bubble and into the 
mixing layer, then the second expansion law holds. 

studied wind-ISM and wind-wind interactions with multi- 
dimensional simulations, showing that some shocks become effectively isothermal, form- 
ing thin and unstable layers that break up into clumps. This multi-dimensional effect is 
important for the expansion rate and radiative emission of the nebulae, and can explain 
some of the observational properties of Wolf-Rayet nebulae. This line of research was 
extended by with the inclusion of photoionizing radiative transfer, 
modelling the H 0 region and wind-bubble evolution for the full stellar lifetime, expanding 
into a uniform ISM. They showed that dissipation and instabilities in expanding layers 
can significantly affect X-ray and UV emission from wind bubbles, potentially resolving 
the disagreement between large predicted X-ray luminosities, and relatively weak X-ray 
emission from observed nebulae. used different stellar evolution 
sequences and higher-resolution simulations, including calculations with and without 
thermal conduction, investigating the structure and X-ray emission from massive-star 
nebulae. They concluded that both thermal conduction and dynamical mixing by hydro- 
dynamic instabilities are affecting the X-ray emission. modelled the 
evolution of the CSM around a massive star including effects of both photoionization 
and stellar wind in 3D simulations with a similar setup. The latest work using as initial 
conditions wind expanding from a star at rest with respect to a uniform ISM is by [Meyer] 
(2020), who used high-resolution 2D simulations to study the CSM of a 60 Mo star 
through various evolutionary phases followed by explosion as a supernova. 

These simulations are very useful for elucidating the physical processes at work in wind 
bubbles, but they assume certain symmetries that are not always present. In particular 


stellar motion or bulk flows in the ISM lead to distorted wind bubbles (Weaver et al.]1977 
2015) that produce a bow shock when the star is moving supersonically 
with respect to the ISM (e.g./van Buren & McCray|1988). Furthermore the turbulent ISM 


is neither static nor homogeneous, and these density and velocity fluctuations also affect 


the shape of a wind bubble (see Fig. |1| taken from {Geen et al. 2021). The combination 


of stellar winds, radiative feedback, stellar clustering, and supernovae leads to a very 
complex multiphase ISM with strong pressure gradients on many scales (Rathjen et al. 


021). 


2. Runaway Stars and bow shocks 
Not long after the first calculations showing that the interaction between Solar wind 


(Baranov et al.]1971) or stellar wind (Weaver et al.)1977) and local ISM could produce a 


bow shock, the bow shock around the closest O star to Earth, ¢ Ophiuchi, was detected in 


nebular emission lines by|Gull & Sofia) (1979). Bow shocks provide an excellent laboratory 


for studying the wind-ISM interaction for a number of reasons: (i) the star has generally 
moved far from its place of birth and there are no other massive stars in the immediate 
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Figure 1. Slice through a 3D simulation of a turbulent molecular cloud, 0.4 Myr after a massive 
star has formed, showing gas density (left) and temperature (right). The hot, low-density and 


asymmetric wind bubble surrounds the star. Reproduced from fig. 4 of “The geometry and 
dynamical role of stellar wind bubbles in photoionized H I regions” Lesser (2021), MNRAS, 
501, 1352. 


vicinity that could induce wind-wind interactions; (ii) the star is typically moving in the 
diffuse (warm) phase of the ISM through a less structured medium than in a molecular 
cloud; and (iii) ram pressure provided by stellar motion compresses the bow shock into 
a structure that is overdense (and more easily observable) and that has a relatively short 
dynamical timescale. 

The structure of a bow shock is well-described in fig. 1 of (1998), 
again consisting of two shocks separated by a contact discontinuity. Stellar motion at 
velocity v, with respect to the ISM leads to an asymmetric ISM ram pressure of pov? in 
the reference frame of the star (po is the ISM gas density). Pressure balance gives the 


characteristic size of the bow shock (Baranov et al.|1971), the standoff distance, Ro, as 


Mu 
Ro = 21 
S don Ju +c?) e 


where M and væ are the mass-loss rate and terminal wind velocity of the star, respec- 
tively. The Alfvén (va) and sound (c) speeds represent the magnetic and thermal pressure 
contributions to the total pressure, but usually the ram pressure is the dominant term. 
The dynamical timescale of the bow shock is tayn = Ro/vs, typically ~ 10° yr for bow 
shocks with Ro < 1pc and v, ~ 30kms~!, much shorter than the lifetime of a massive 
star. 

This equation provides a measurement of M that is independent from stellar-atmosphere 


spectral-line studies (e.g. (Gull & Sofia||1979), if the other quantities on the right-hand 
side of the equation are well constrained. |Gvaramadze et al.| (2012) were able to use 


the properties of the H II region around the runaway O star ¢ Oph to constrain pp and 
thereby measure M x 2.2 x 1078 Mo yr™!1, below the estimate derived from optical lines 
in the atmosphere and the (2000) theoretical estimate, but well above the 
estimate derived from UV lines [Marcolino et al.|{2009). use 
Eq. [2.1] to measure M for a sample of bow shocks around 20 O stars by making the 
assumptions that all outer shocks are adiabatic with a density jump of 4x, that po can 
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be accurately measured from far-IR dust emission, and all stars are moving through the 


ISM with v, = 30kms~!. This was followed up by |Kobulnicky et al.| (2019) who used 


Gaia DR2 data to measure v, for a larger sample of stars and thereby obtaining more 


accurate results. made a careful analysis of the statistical and 
systematic uncertainties of measuring M from bow-shock observations, also comparing 
their method with that of (2018). They find that M measurements 
for individual sources may have large uncertainties (from uncertain dust properties and 


shocked ISM pressure support), but statistically for a large group of sources the methods 
are promising. 


2.1. Multiwavelength emission from bow shocks 


The first detection of a significant number of bow shocks was made with the IRAS 
observatory by in the mid-IR (see also [van Buren et al.| 
(1995). The broadband IR emission is thermal radiation from interstellar dust grains 
photo-heated by the intense UV-radiation field of the massive star that drives the bow 
shock. The AKARI observatory obtained higher-resolution IR images of some bow shocks 


(Ueta et al.[2008), but a major advance came with the Spitzer and WISE missions. The 
E-BOSS "E (2012| contains 73 confirmed and candidate bow 
shocks, and (2016) compiled a larger catalog of 709 candidate bow 
shocks selected by mid-IR. morphology. (2014) showed that bow shocks are 
most luminous in the mid-IR, (see also [Henney & Arthurl/2019), 
explaining why these surveys have been so successful. Asymmetric wind bubbles with 
subsonic relative motion between star and ISM can also produce bright mid-IR arcs 
(Mackey et al.[2016). 

Despite the first detection of a bow shock being in optical lines, this has proven to 
be a difficult method for detection because the massive stars are so optically bright that 
identifying nearby nebular emission is challenging (see [Gull & Sofial[1979). A few other 
detections have been made: HD 165319 appears to show a bow shock in narrow-band 
Ha imaging (Gvaramadze & Bomans]2008), and the red supergiant IRC -10414 also has 
a bow shock detected in Ha imaging (Gvaramadze et al./2014). It seems likely that the 
bow shocks of Betelgeuse (Mohamed et al.||2012) and € Ophiuchi 


could be detectable in nebular optical lines if the exceptionally bright stellar emission 
could be masked. 

The first radio detection of a bow shock was by (2010), around the run- 
away O supergiant BD+43 3654, further studied by (2021) and {Moutzouril 
over a wider frequency range. There is significant interest in radio studies 
at present because of the rapid advances in instrumentation and also the potential to 
detect both thermal (Bremsstrahlung) and non-thermal (synchrotron) radiation. From 
its radio spectrum it seems that the bow shock of BD+43 3654 exhibits both. 
also detected the Bubble Nebula, NGC 7635, as predicted NEE 
(2019), and the detections of this and BD+43 3654 are shown in Fig.|3} Some low signif- 
icance radio emission was found in NVSS survey data for E-BOSS bow shocks by 
(2015), and some of these are now confirmed radio emitters with ASKAP observa- 
tions (Van den Eijnden et al. (20221. The bow shock of Vela X-1 was also detected with 
MeerKAT (van den Eijnden et al.|/|2022). Given these recent discoveries and the promise 
of wide-field, broadband surveys with the new interferometric arrays, I anticipate many 
more radio detections of bow shocks in the coming years. 

Thermal bremsstrahlung gives a direct measure of the emission measure, allowing to 
constrain the gas density in both the pre-shock and shocked ISM in the bow shock. 
Obtaining gas density from IR emission has significant uncertainties arising from dust 
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Figure 2. Four bow shocks from the catalog o Dee 2015), imaged by WISE. Red: band 
4, 22.2 um. Green: band 3, 12.1 um. Blue: band 4 um. redit: Peri, Benaglia & Isequilla, 
2015, Astronomy & Astrophysics, 578, A45, fig. 5 


composition and grain size distribution (Pavlyuchenkov et al./2013}|Mackey et al. (20101. 


and nebular line emission can have uncertainties from patchy extinction along the line 
of sight. Radio has the distinct advantage that is it unaffected by extinction. 


With a detection of synchrotron radiation at radio wavelengths one can make model- 
dependent predictions for synchrotron and inverse-Compton emission at X-ray and y-ray 


energies (del Valle & Romero 2012 del Palacio et al.|/2018; 2022). 
Searches for non-thermal X-rays (Toala et al.!2016||2017) and y-rays (Schulz et al./2014 
H.E.S.S. Collaboration et al.|/2018) from bow shocks have so far only obtained upper 


limits, although this is not altogether surprising given the predictions and sensititivies 
of current instruments. Detection of high-energy radiation may need to wait for the next 


generation of instruments for bow shocks around single stars (Moutzouri et al.|/2022). 


6 Jonathan Mackey 


T T T T 0.0040 T T T T 0.040 
gg SEIST I 0.035 
iio J$ 0.003: 0.035 
0.0030 0.030 
Lt | 
s s 0.0025 0.025 
3 Dt é W 4 
S "SE, 0.0020 INS 4F} 0.020 
Q a e Te af "` 
R ji S o 
A "oi l K : J 10.0015 0.015 
LA : 4 
10'b i J 
© © |f 40.0010 i 0.010 
CG 
43°55'b — “Jr 70.0005 S eg 0.005 
` IEN "4 
! : © > 0.0000 ` ! ! S 0.000 
2034™305 005 331305 00s ` 2321™205 ` or 20™405 205 i 
RA (J2000) RA (J2000) 


Figure 3. Radio detection of the bow shocks around BD+43 3654 (left) and BD+60 2522 (right), 
obtained with the VLA at 4-8 GHz (greyscale and grey contours). Blue contour shows WISE 
band 4 emission in mid-IR, and cyan contours show 1.4GHz emission from the NVSS survey. 
Credit: Moutzouri et al., 2022, Astronomy & Astrophysics, 663, A80, fig. 1. 


3. The wind-ISM boundary layer — theoretical models 


The contact discontinuity between shocked stellar wind and interstellar material is 
very strong for hot stars with fast winds. Typical post-shock temperature of the shocked 
wind is Tpg ~ 10° K for an adiabatic shock, easily obtainted from the Rankine-Hugoniot 
jump conditions: 


2 
Ush 
Tps =~ 5.6 x 107 K | = — 4 3.1 
SN Lal GE 


where an exact expression depends on the gas composition and ionization state. Given 
a photoionized ISM temperature of Tig), œ 8000 K, this implies a temperature jump of 
a factor of 104 across the contact discontinuity, and an associated density jump of the 
same factor in the opposite sense. For comparison this is similar to the density ratio of 
air and rock. It is easy to see that mixing a small volume of interstellar gas across the 
contact discontinuity can have a large effect on the density and temperature (and hence 
radiative emissivity at different energies) of the wind bubble or bow shock. 

Electronic thermal conduction transports heat from the hot and 
low-density gas phase across a contact discontinuity into a boundary layer of cooler and 
denser gas. The thickness of this boundary layer is determined by the mean free path 
of electrons, or the gyroradius in the presence of a magnetic field component parallel 
to the discontinuity. The heated boundary layer expands, resulting in a region with 
strong temperature and density gradients rather than a well-defined contact discontinuity. 


This can be seen in hydrodynamic bow shock simulations by (1998) 
and (2014), where the structure of the bow shock is signficantly modified 


by thermal conduction. |Meyer et al.) (2017) showed, however, that the inclusion of an 
interstellar magnetic field dramatically reduces the effects of thermal conduction and, 
given that the ISM is pervaded by magnetic fields, one may deduce that the hydrodynamic 


models of wind bubbles (Weaver et al.||1977) and bow shocks (Comerén & Kaper|/1998 
2014) significantly overestimate the effects of thermal conduction. 


An alternative physical process that can thicken and distort the wind-ISM bound- 
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Figure 4. Gas density (gcm~?, upper half-plane) and temperature (K, lower half-plane) on a 
logarithmic colour scale for high-resolution 2D axisymmetric simulations of bow shock models 
for the Bubble Nebula, NGC 7635. Credit: Green et al., 2019, Astronomy & Astrophysics, 625, 
A4, fig. 2. 


ary layer is often referred to as dynamical or turbulent mixing (e.g. 
applied to supernova remnants). This is generally driven by Kelvin-Helmholtz insta- 
bility (KHI) because the hot phase has a sound speed up to 100x that of the warm 
ISM phase, and so the velocity shear at the interface is typically large. It cannot be 
captured in 1D simulations, but 2D simulations with an unstable wind-ISM interface 
(due to e.g. Rayleigh-Taylor instability or non-linear thin-shell instability) or with an 
anisotropic external pressure E generate such shear flows (Toalá & Arthur (Toalé & Arthur|2011} 
[Mackey et_al.|/2015 2015} Green et al.|/2019). This KHI-induced turbulent mixing is clearly 
shown in Fig. Hl (taken from Seen et al. Green et _al.][2019), where waves initiated near the apex 
of the bow shock grow to non-linear amplitudes and the resulting vortices strongly mix 
wind and ISM material in the wake behind the bow shock. 

If thermal conduction is not explicitly modelled then the degree of mixing is dependent 
on numerical resolution through both numerical diffusivity and the degree to which the 
numerical scheme has the resolution to capture the development of KHI (e.g. 
et al.[2022). Recently, [Lancaster ct al] (2021) and [Lancaster SU 2021b) developed a 
theory of turbulent mixing in wind bubbles expanding into clouds with a fractal density 
structure, arguing that this should lead to strongly cooled wind bubbles. The degree of 
energy dissipation is a key constraint for larger scale simulations of the ISM of galaxies, 
for which it is currently not clear how efficient wind feedback (implemented in a sub-grid 


model) should be (e.g. 2022). 


4. The wind-ISM boundary layer — observations 


One way to constrain these processes observationally is to study the thermal X-ray or 
UV emission from the wind-ISM interface. The brightest X-ray emitting wind bubbles 


are around Wolf-Rayet (WR) stars (Toalá et al 2018). These are wind-wind interactions 


where the fast wind of the WR star sweeps up the slow wind (García-Segura et al.]1996b 
or mass ejection from binary systems (Jiménez-Hernandez et al.||2021). Here there is 
usually considerable uncertainty about the density structure of the wind in the previous 
evolutionary phase. UV lines from intermediate temperature gas at the hot-cold interface 
were also detected by {Boroson et al. within a WR ring nebula and interpreted as 
a conduction front. 
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So far the only clear detection of diffuse X-rays from stellar wind of a single main- 
sequence star is very weak emission from within the bow shock around Ç Oph (Toalé| 
2022). Diffuse X-rays were also predicted (Mackey et al.|2015) 
and later detected (Townsley et al. 20191 from the massive star forming region RCW 120, 
which likely includes a contribution from the O star that is the main driver of the nebula’s 
evolution. A quantitative comparison between theory and observation for this object has 
not been made. 

On larger scales, X-ray emission has been detected from a number of massive star- 


forming regions (Townsley et al./2018). The energy content of the hot phase was derived 


from observational data and has been compared with the energy input from winds by 
(2014), finding that most of the input energy cannot be accounted for. 
This is strong observational evidence for the effectiveness of energy transport across the 
contact discontinuity by e.g. turbulent mixing, although ‘leakage’ of energy into the low- 
density coronal gas surrounding star-forming regions also could not be excluded 


E Pittard[2014). 


5. 3D MHD models of bow shocks and astrospheres 


We have seen that magnetic fields are important for moderating the effects of thermal 
conduction at the contact discontinuity of bow shocks, and it is also known that the 
operation of KHI is markedly different in the presence of a magnetic field 
[Keppens et al.|[1999). Axisymmetric simulations in 2D do not allow a general 
magnetic field configuration because only ISM fields that are parallel with (or anti- 
parallel to) the velocity vector of the star are permitted. There is therefore a need for 3D 
MHD simulation of bow shocks, and a number of groups have worked on this over the 
past few years, developing efficient methods that enable 3D simulations with reasonable 
computational resources. 

The computationally cheapest calculations are for stars with slow winds (because the 
timestep is inversely proportional to the fastest speed on the domain), and so the first 
3D hydrodynamic simulations were for isothermal calculations of winds from cool stars 


(Blondin & Koerwer}1998). These showed very strong dynamical instability that quickly 
grew to nonlinear amplitudes. More detailed calculations by |Mohamed et al.) (2012) 


applied to the red supergiant Betelgeuse explored the stability of the bow shock as 
a function of space velocity of the star, with a detailed radiative cooling prescription. 
Recently [Meyer et a) added magnetic fields to these models of bow shocks around 
cool stars, showing that the field has a strong stabilising influence on the bow shock. 

Simulations of bow shocks around O stars are much more challenging because of the 
larger Mach number of the shocks (affecting stability of the scheme) and the timestep 
constraint (increasing the computational cost). [Scherer et al (2020) developed a 3D MHD 
scheme using a spherical coordinate system that naturally has varying spatial resolution 
and hence is quite efficient. They studied the differences between hydrodynamic and MHD 
simulations for a bow shock similar to that around the massive stars \ Cep. 
used this model to study the effects of ISM inhomogeneities on the structure 
and observable properties of bow shocks. 


Using a Cartesian coordinate system with static mesh-refinement,|Mackey et al.| (2021) 


introduced 3D MHD simulations of bow shocks with the software PION. This was used by 
(2022), who presented the first 3D MHD simulations dedicated to modelling 
the bow shock of ¢ Oph. This spectacular system (see Fig. provides arguably the 
cleanest laboratory for testing mixing processes at the wind-ISM boundary. It is a very 
nearby massive star (136 pc) with well-characterised proper motion, moving through a 
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Figure 5. Left: The bow shock of ¢ Ophiuchi, seen in infrared (red and green, from 
Spitzer Space Telescope) and X-rays (blue, Chandra X-ray Observatory). Image credit: X-ray: 
NASA/CXC/Dublin Inst. for Advanced Studies/S. Green et al.; Infrared: NASA/JPL/Spitzer. 
Right: synthetic X-ray emission from 3D MHD simulation of the bow shock, from {Green et al.] 
(Astronomy & Astrophysics, 665, A35, fig. 13). 


large photoionized H 1 region and driving a bright and well-resolved bow shock. Cru- 
cially, diffuse X-ray emission has been detected with Chandra (Toalé et al.//2016). This 
source is also far out of the Galactic Plane, which can be valuable for multi-wavelength 
studies because background emission is less of an issue. Becuase of rapid stellar rota- 
tion the radial velocity of the star is actually poorly constrained, and this introduced 
some uncertainty to the modelling. Nevertheless the authors were able to show that the 
X-ray emission predicted by the simulated bow shock is somewhat weaker than the ob- 
served emission. There are uncertainties, particularly the role of instabilities and effects 
of limited numerical resolution, but this comparison of 3D MHD simulations with high- 
resolution observations is a promising avenue for obtaining detailed information on the 
physical processes operating at shocks and contact discontinuities in the hot and warm 
ISM phases. 


6. Outlook 


The interaction of stellar winds with the ISM is important for interpreting observations 
of circumstellar nebulae, but also has broader application to a number of other areas of 
astrophysics. The effectiveness of thermal conduction in astrophysical plasmas is so far 
only experimentally measured in the Solar wind (e.g.|Bale et al 12012. and the degree to 
which it is inhibited by both large-scale coherent and small-scale turbulent magnetic fields 
can potentially be investigated at the wind-ISM interface. It is important to constrain 
mixing of interstellar matter into stellar-wind bubbles and how this affects the efficiency 
of mechanical feedback from winds to the ISM, with consequences for models of galaxy 
formation and evolution. The wind-ISM interaction sets the initial conditions into which a 
supernova blast wave expands, which can have important consequences for the evolution 


of supernova remnants (Das et al 12022. Bow shocks and wind bubbles can be a useful 
laboratory for studying particle acceleration in shocks (Benaglia et al.||2010} |H.E.S.S. 
esses (del 


Collaboration et al 20221 and associated particle transport and radiation proc 
Valle & Pohl/2018). 


This review has focussed on wind bubbles and bow shocks around single stars, but 
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I want to emphasise that the coming decade will see great advances in modelling and 
observing wind-ISM and wind-wind interactions of binary and higher-order multiple star 
systems. Unlike low-mass stars, almost all massive stars begin their lives in binaries 
and the majority will undergo interaction with a companion during their lifetime 
fet al DO, Colliding-wind binaries are one of the few astrophysical systems where 


time-dependent shock physics can be probed (Pittard & Dougherty|2006}|H.E.S.S. Col- 
2020), and they are bright Galactic sources across the electromagnetic 


spectrum up to TeV gamma rays. 

Progress in our understanding of wind-ISM interaction depends on both new data 
and more detailed models. In this regard we can anticipate breakthroughs within the 
next decade. The rapid improvement in large field-of-view radio interferometry driven 
by SKA pathfinder instruments is leading to detection of both thermal and non-thermal 
radio emission from circumstellar nebulae and bow shocks. This gives important in- 
sights into gas density and thermal state, and the population of non-thermal particles. 
In the next few years the Cherenkov Telescope Array (CTA) will come online, provid- 
ing unprecedented sensitivity to TeV gamma-ray emission, with potential detections of 
populations of very-high-energy particles accelerated in wind bubbles and bow shocks, 
and detection of a large population of Galactic binary systems, potentially including the 
colliding-wind binaries. Looking further ahead, the ATHENA mission promises huge sen- 
sitivity improvements compared with current X-ray telescopes, leading to much better 
characterisation of the hot thermal plasma of the shocked stellar-wind, and potentially 
detection of non-thermal emission via synchrotron or inverse-Compton radiation. 

At the same time, rapid improvements in software for astrophysical fluid dynamics, 
including open-source community projects, mean that high-fidelity simulations of bow 
shocks and wind bubbles can now be performed. In comparison with new high-resolution 
datasets coming from observations across the electromagnetic spectrum, these models 
have significant power to constrain uncertain physical parameters, discriminate between 
different physical models, and in general give a deeper and clearer understanding of 
astrophysical gas dynamics. 
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